function[y] = fun(params,varargin)
if nargin>1
    saveas = char(varargin(1));
else
    saveas = 'Baseline';
end
%fprintf(saveas);
%fprintf('\n');

psi       = 0.03;
alph      = 0.25; %0.408; %M_data(12,1);   % the experience advantage in getting listings
bet       = 0.9;   % Discount factor
M_data = xlsread('data_for_model/save_calibrationcontall_experience_l1', 'bust'); 
s(1) = M_data(11,1);
M_data = xlsread('data_for_model/save_calibrationcontall_experience_l1', 'medium'); 
s(2) = M_data(11,1);
M_data = xlsread('data_for_model/save_calibrationcontall_experience_l1', 'boom'); 
s(3) = M_data(11,1);
S_state   = [s(1) s(2) s(3)];  % number of sellers (bust/medium/boom)
kh        = 3;     % aggregate states
ka        = 1;     % bins for ability
ke        = 50;    % bins for experience

l_bar   = 60;    % capacity constraint, will only matter if it's above the experience level? hmm...


%% Set up of states
kw    = kh^2;     % number of histories of aggregate states (1-9)
e_bar   = min(l_bar,ke); % maximum level of experience that matters
State = combvec(1:1:ke,(1:1:ka)./ka,1:1:kw); % ability vector translates into the multiplier
State(4,:) = mod(State(3,:),kh)+kh.*(mod(State(3,:),kh)==0); % actual aggregate state (1,2,3)
ks = max(size(State)); % state space
init = 1/ka*(State(1,:)==1); % initial distr. uniform across ability DOES THIS MATTER? (SG 2023)

%% State probability function
% expand the State matrix
S_state     = repmat(S_state,1,kw/kh); % number of sellers in eah state history

%% one way: assume one state markov
% the following procedure works for one lag only!
load data_for_model/transition_matrices T TT
T_h_vec = kh.*(mod(1:1:kw,kh)==0);
h_vec = kh.*(mod(State(3,:),kh)==0);
T_h_vec_lag = kh.*(mod(ceil((1:1:kw)./kh),kh)==0);
for i=1:kh-1 
  h_vec   = h_vec+i.*(mod(State(3,:),kh)==i);
  T_h_vec = T_h_vec+i.*(mod(1:1:kw,kh)==i);
  T_h_vec_lag =  T_h_vec_lag+i.*(mod(ceil((1:1:kw)./kh),kh)==i);
end
T_small = T;
ind = sub2ind(size(T_small),repmat(T_h_vec',1,kw),repmat(T_h_vec,kw,1));
T = reshape(T_small(ind),kw,kw).*(repmat(T_h_vec',1,kw)==repmat(T_h_vec_lag,kw,1));

% compute the steady state aggregate distribution
[V,D] = eig(T'); % Find eigenvalues and left eigenvectors of A
[~,ix] = min(abs(diag(D)-1)); % Locate an eigenvalue which equals 1
q = V(:,ix)'; % The corresponding row of V' will be a solution
q = q/sum(q); % Adjust it to have a sum of 1
T(q==0,:)=0;
q_vec = q(State(3,:));
license_limit = 1000000000000000;

y=model_passive_learn_dyn_prices_poisson(params, ke, State, kh, ka, T, ks, kw, l_bar, init, psi, alph, bet, S_state, q_vec, q,license_limit,saveas,0,0);

